Historical Contingency of Host-Microbiome Response to Perturbation

Overview

Setup

(hide)

Click on tabs to display additional information.

Libraries

# Set seed for reproducibility
set.seed(42)

# Load required libraries for DESeq2 analysis and visualization
library(DESeq2)
library(SummarizedExperiment)
library(ggplot2)
library(pheatmap)
library(RColorBrewer)
library(ggrepel)
library(tidyverse)
library(VennDiagram)
library(gridExtra)
library(BiocParallel)
library(gt)
library(broom)
library(car)
library(emmeans)
library(multcomp)
library(MASS)
library(rstatix)
library(coin)

Plotting

# Define treatment order and color palette for plotting (original naming convention)
treatment_order_plot <- c(
  "A- T- P-",  # Control
  "A- T- P+",  # Parasite
  "A+ T- P-",  # Antibiotics
  "A+ T- P+",  # Antibiotics_Parasite
  "A- T+ P-",  # Temperature
  "A- T+ P+",  # Temperature_Parasite
  "A+ T+ P-",  # Antibiotics_Temperature
  "A+ T+ P+"   # Antibiotics_Temperature_Parasite
)

# Custom color palette matching treatment order
treatment_colors <- c(
  "#1B9E77",  # A- T- P- (Control)
  "#D95F02",  # A- T- P+ (Parasite)
  "#7570B3",  # A+ T- P- (Antibiotics)
  "#E7298A",  # A+ T- P+ (Antibiotics_Parasite)
  "#66A61E",  # A- T+ P- (Temperature)
  "#E6AB02",  # A- T+ P+ (Temperature_Parasite)
  "#A6761D",  # A+ T+ P- (Antibiotics_Temperature)
  "#666666"   # A+ T+ P+ (Antibiotics_Temperature_Parasite)
)

# Create named vector for color scale
treatment_color_scale <- setNames(treatment_colors, treatment_order_plot)

# Define treatment order for DESeq2 (keep original DESeq2-compatible names)
treatment_order_deseq <- c(
  "A_minus_T_minus_P_minus",  # Control
  "A_minus_T_minus_P_plus",   # Parasite
  "A_plus_T_minus_P_minus",   # Antibiotics
  "A_plus_T_minus_P_plus",    # Antibiotics_Parasite
  "A_minus_T_plus_P_minus",   # Temperature
  "A_minus_T_plus_P_plus",    # Temperature_Parasite
  "A_plus_T_plus_P_minus",    # Antibiotics_Temperature
  "A_plus_T_plus_P_plus"      # Antibiotics_Temperature_Parasite
)

# Create mapping between DESeq2 names and plotting names
treatment_mapping <- setNames(treatment_order_plot, treatment_order_deseq)

# Mapping from DESeq2 contrast names to plotting-friendly names
contrast_to_plotname <- c(
  "A_minus_T_minus_P_minus_vs_A_minus_T_minus_P_minus" = "A- T- P-",
  "A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus" = "A- T- P+",
  "A_plus_T_minus_P_minus_vs_A_minus_T_minus_P_minus" = "A+ T- P-",
  "A_plus_T_minus_P_plus_vs_A_minus_T_minus_P_minus" = "A+ T- P+",
  "A_minus_T_plus_P_minus_vs_A_minus_T_minus_P_minus" = "A- T+ P-",
  "A_minus_T_plus_P_plus_vs_A_minus_T_minus_P_minus" = "A- T+ P+",
  "A_plus_T_plus_P_minus_vs_A_minus_T_minus_P_minus" = "A+ T+ P-",
  "A_plus_T_plus_P_plus_vs_A_minus_T_minus_P_minus" = "A+ T+ P+"
)

Functions

# Function to convert DESeq2 treatment names to plotting names
convert_treatment_names <- function(deseq_names) {
  return(treatment_mapping[deseq_names])
}

# Function to monitor system resources during analysis
monitor_performance <- function() {
  # Get memory usage
  mem_usage <- gc()
  
  # Get CPU usage (if possible)
  cpu_info <- tryCatch({
    if (Sys.info()["sysname"] == "Darwin") {
      # macOS
      cpu_cmd <- "top -l 1 -n 0 | grep 'CPU usage'"
      cpu_output <- system(cpu_cmd, intern = TRUE)
      cpu_output
    } else if (Sys.info()["sysname"] == "Linux") {
      # Linux
      cpu_cmd <- "top -bn1 | grep 'Cpu(s)'"
      cpu_output <- system(cpu_cmd, intern = TRUE)
      cpu_output
    } else {
      "CPU monitoring not available on this system"
    }
  }, error = function(e) {
    "CPU monitoring not available"
  })
  
  cat("=== Performance Monitor ===\n")
  cat("Memory usage:\n")
  print(mem_usage)
  cat("CPU info:\n")
  cat(cpu_info, "\n")
  cat("Active parallel workers:", bp_param$workers, "\n")
  cat("========================\n")
}

# Function to optimize memory usage
optimize_memory <- function() {
  cat("Optimizing memory usage...\n")
  
  # Force garbage collection
  gc(verbose = FALSE)
  
  # Clear unnecessary objects
  if (exists("temp_objects")) {
    rm(temp_objects)
  }
  
  cat("Memory optimization complete\n")
}

# Function to save DESeq2 results
save_deseq_results <- function(results_list, analysis_name) {
  # Create results directory if it doesn't exist
  results_dir <- here::here("Code", "Analysis", "DiffExpGene", "Saved_Results")
  dir.create(results_dir, recursive = TRUE, showWarnings = FALSE)
  
  # Save the results list
  results_file <- file.path(results_dir, paste0(analysis_name, "_deseq_results.rds"))
  saveRDS(results_list, results_file)
  
  cat("Results saved to:", results_file, "\n")
}

# Function to load DESeq2 results
load_deseq_results <- function(analysis_name) {
  results_dir <- here::here("Code", "Analysis", "DiffExpGene", "Saved_Results")
  results_file <- file.path(results_dir, paste0(analysis_name, "_deseq_results.rds"))
  
  if (file.exists(results_file)) {
    cat("Loading existing results for:", analysis_name, "\n")
    return(readRDS(results_file))
  } else {
    cat("No existing results found for:", analysis_name, "\n")
    return(NULL)
  }
}

# Function to check if results exist
results_exist <- function(analysis_name) {
  results_dir <- here::here("Code", "Analysis", "DiffExpGene", "Saved_Results")
  results_file <- file.path(results_dir, paste0(analysis_name, "_deseq_results.rds"))
  return(file.exists(results_file))
}

# Function to clean up saved results
cleanup_saved_results <- function(analysis_names = NULL) {
  results_dir <- here::here("Code", "Analysis", "DiffExpGene", "Saved_Results")
  
  if (is.null(analysis_names)) {
    # Remove all saved results
    if (dir.exists(results_dir)) {
      unlink(results_dir, recursive = TRUE)
      cat("All saved results removed\n")
    }
  } else {
    # Remove specific analysis results
    for (analysis_name in analysis_names) {
      results_file <- file.path(results_dir, paste0(analysis_name, "_deseq_results.rds"))
      if (file.exists(results_file)) {
        file.remove(results_file)
        cat("Removed saved results for:", analysis_name, "\n")
      }
    }
  }
}

# Function to run DESeq2 analysis for a single comparison
run_deseq_comparison <- function(dds, contrast_name, contrast_vector) {
  cat("Running analysis for:", contrast_name, "\n")
  
  # Set seed for reproducibility within parallel processes
  set.seed(42)
  
  # Optimize DESeq2 settings for speed
  # Use faster estimation method for large datasets
  if (nrow(dds) > 10000) {
    # For large datasets, use faster estimation
    dds_comparison <- DESeq(dds, 
                           parallel = TRUE,
                           fitType = "local",  # Faster than "parametric"
                           betaPrior = FALSE,  # Disable beta prior for speed
                           quiet = FALSE)
  } else {
    # For smaller datasets, use standard settings
    dds_comparison <- DESeq(dds, 
                           parallel = TRUE,
                           quiet = FALSE)
  }
  
  # Get results with optimized settings
  res_comparison <- results(dds_comparison, 
                           contrast = contrast_vector,
                           parallel = TRUE,
                           alpha = 0.05)  # Set significance threshold
  
  # Add gene names if available
  if ("gene_name" %in% names(mcols(dds_comparison))) {
    res_comparison$gene_name <- mcols(dds_comparison)$gene_name
  }
  
  # Clean up memory
  rm(dds_comparison)
  gc()
  
  return(res_comparison)
}

# Function to run all comparisons for a specific analysis with enhanced parallelization
run_analysis_set <- function(analysis_name, analysis_config, force_recompute = FALSE) {
  cat("\n=== Running", analysis_name, "===\n")
  cat("Description:", analysis_config$description, "\n")
  
  # Check if results already exist and force_recompute is FALSE
  if (!force_recompute && results_exist(analysis_name)) {
    cat("Loading existing results for", analysis_name, "\n")
    return(load_deseq_results(analysis_name))
  }
  
  # Create a copy of the dataset for this analysis
  dds_analysis <- dds_filtered
  
  # Update design if needed
  design(dds_analysis) <- analysis_config$design
  
  # Pre-estimate size factors and dispersions for better parallelization
  cat("Pre-estimating size factors and dispersions...\n")
  dds_analysis <- estimateSizeFactors(dds_analysis)
  dds_analysis <- estimateDispersions(dds_analysis, parallel = TRUE)
  
  # Run all contrasts for this analysis
  results_list <- list()
  
  # Use parallel processing for multiple contrasts
  if (length(analysis_config$contrasts) > 1) {
    cat("Running", length(analysis_config$contrasts), "contrasts in parallel...\n")
    
    # Create a function for parallel execution
    run_single_contrast <- function(contrast_info) {
      contrast_name <- contrast_info$name
      contrast_vector <- contrast_info$vector
      
      # Check if both treatment levels exist in the data
      if (all(contrast_vector[2:3] %in% levels(colData(dds_analysis)$Treatment))) {
        return(list(name = contrast_name, result = run_deseq_comparison(dds_analysis, contrast_name, contrast_vector)))
      } else {
        cat("Warning: Contrast", contrast_name, "skipped - treatment levels not found in data\n")
        return(list(name = contrast_name, result = NULL))
      }
    }
    
    # Prepare contrast information for parallel processing
    contrast_info_list <- lapply(names(analysis_config$contrasts), function(name) {
      list(name = name, vector = analysis_config$contrasts[[name]])
    })
    
    # Run contrasts in parallel
    parallel_results <- BiocParallel::bplapply(contrast_info_list, run_single_contrast)
    
    # Organize results
    for (result in parallel_results) {
      if (!is.null(result$result)) {
        results_list[[result$name]] <- result$result
      }
    }
    
  } else {
    # For single contrast, run directly
    for (contrast_name in names(analysis_config$contrasts)) {
      contrast_vector <- analysis_config$contrasts[[contrast_name]]
      
      # Check if both treatment levels exist in the data
      if (all(contrast_vector[2:3] %in% levels(colData(dds_analysis)$Treatment))) {
        results_list[[contrast_name]] <- run_deseq_comparison(dds_analysis, contrast_name, contrast_vector)
      } else {
        cat("Warning: Contrast", contrast_name, "skipped - treatment levels not found in data\n")
      }
    }
  }
  
  # Save the results
  save_deseq_results(results_list, analysis_name)
  
  # Clean up memory
  rm(dds_analysis)
  gc()
  
  return(results_list)
}

# Function to format p-values for gt tables
fmt_pvalues <- function(gt_tbl, columns, threshold = 0.001, decimals = 3) {
  gt_tbl |>
    gt::fmt_number(columns = all_of(columns), decimals = decimals) |>
    gt::fmt(
      columns = all_of(columns),
      rows = Reduce(`|`, lapply(columns, function(col) gt_tbl$`_data`[[col]] < threshold)),
      fns = function(x) rep(paste0("<", threshold), length(x))
    )
}

# Function to save results for each analysis
save_analysis_results <- function(analysis_name, results_list) {
  cat("Saving results for:", analysis_name, "\n")
  
  # Create subdirectory for this analysis
  analysis_dir <- file.path(output_dir, analysis_name)
  dir.create(analysis_dir, recursive = TRUE, showWarnings = FALSE)
  
  # Save each contrast result
  for (contrast_name in names(results_list)) {
    res <- results_list[[contrast_name]]
    
    # Convert to data frame
    res_df <- as.data.frame(res)
    
    # Add gene names if available
    if ("gene_name" %in% names(res_df)) {
      res_df <- res_df %>% 
        dplyr::select(gene_name, everything()) %>%
        dplyr::rename(gene = gene_name)
    }
    
    # Save full results
    full_filename <- file.path(analysis_dir, paste0(contrast_name, "_full_results.tsv"))
    write_tsv(res_df, full_filename)
    
    # Save significant results (padj < 0.05)
    sig_df <- res_df %>% 
      dplyr::filter(padj < 0.05) %>%
      dplyr::arrange(padj)
    
    if (nrow(sig_df) > 0) {
      sig_filename <- file.path(analysis_dir, paste0(contrast_name, "_significant_results.tsv"))
      write_tsv(sig_df, sig_filename)
      
      # Create summary table
      summary_table <- data.frame(
        Contrast = contrast_name,
        Total_Genes = nrow(res_df),
        Significant_Genes = nrow(sig_df),
        Upregulated = sum(sig_df$log2FoldChange > 0, na.rm = TRUE),
        Downregulated = sum(sig_df$log2FoldChange < 0, na.rm = TRUE),
        stringsAsFactors = FALSE
      )
      
      # Save summary
      summary_filename <- file.path(analysis_dir, paste0(contrast_name, "_summary.tsv"))
      write_tsv(summary_table, summary_filename)
    }
  }
  
  # Create combined summary for the analysis
  all_summaries <- list()
  for (contrast_name in names(results_list)) {
    res <- results_list[[contrast_name]]
    res_df <- as.data.frame(res)
    sig_df <- res_df %>% dplyr::filter(padj < 0.05)
    
    all_summaries[[contrast_name]] <- data.frame(
      Contrast = contrast_name,
      Total_Genes = nrow(res_df),
      Significant_Genes = nrow(sig_df),
      Upregulated = sum(sig_df$log2FoldChange > 0, na.rm = TRUE),
      Downregulated = sum(sig_df$log2FoldChange < 0, na.rm = TRUE),
      stringsAsFactors = FALSE
    )
  }
  
  combined_summary <- bind_rows(all_summaries)
  combined_filename <- file.path(analysis_dir, paste0(analysis_name, "_combined_summary.tsv"))
  write_tsv(combined_summary, combined_filename)
  
  # Display summary table
  cat("\nSummary for", analysis_name, ":\n")
  print(combined_summary %>% gt::gt())
}

Import Data

# Read in metadata sheet
Metadata <- readxl::read_excel(here::here("Data", "Transcriptomics", "Data", "ROL_MajorExperiment__MetadataSheet__Corrected_05092025.xlsx"))

# Read in SummarizedExperiment object
se <- readRDS(here::here("Data", "Transcriptomics", "Results", "rnaseq", "051525", "salmon", "salmon.merged.gene_counts_length_scaled.SummarizedExperiment.rds"))

# Print basic information about the dataset
cat("Number of genes:", nrow(se), "\n")

## Number of genes: 39447

cat("Number of samples:", ncol(se), "\n")

## Number of samples: 72

cat("Sample names:", colnames(se), "\n")

## Sample names: TS047_RoL_RNA_110 TS047_RoL_RNA_112 TS047_RoL_RNA_114 TS047_RoL_RNA_116 TS047_RoL_RNA_119 TS047_RoL_RNA_136 TS047_RoL_RNA_142 TS047_RoL_RNA_145 TS047_RoL_RNA_146 TS047_RoL_RNA_16 TS047_RoL_RNA_170 TS047_RoL_RNA_175 TS047_RoL_RNA_177 TS047_RoL_RNA_20 TS047_RoL_RNA_200 TS047_RoL_RNA_204 TS047_RoL_RNA_229 TS047_RoL_RNA_231 TS047_RoL_RNA_258 TS047_RoL_RNA_262 TS047_RoL_RNA_286 TS047_RoL_RNA_288 TS047_RoL_RNA_289 TS047_RoL_RNA_291 TS047_RoL_RNA_317 TS047_RoL_RNA_318 TS047_RoL_RNA_319 TS047_RoL_RNA_322 TS047_RoL_RNA_326 TS047_RoL_RNA_328 TS047_RoL_RNA_347 TS047_RoL_RNA_355 TS047_RoL_RNA_377 TS047_RoL_RNA_381 TS047_RoL_RNA_407 TS047_RoL_RNA_412 TS047_RoL_RNA_439 TS047_RoL_RNA_442 TS047_RoL_RNA_46 TS047_RoL_RNA_468 TS047_RoL_RNA_469 TS047_RoL_RNA_477 TS047_RoL_RNA_478 TS047_RoL_RNA_498 TS047_RoL_RNA_500 TS047_RoL_RNA_502 TS047_RoL_RNA_506 TS047_RoL_RNA_51 TS047_RoL_RNA_526 TS047_RoL_RNA_528 TS047_RoL_RNA_530 TS047_RoL_RNA_532 TS047_RoL_RNA_559 TS047_RoL_RNA_561 TS047_RoL_RNA_562 TS047_RoL_RNA_586 TS047_RoL_RNA_616 TS047_RoL_RNA_618 TS047_RoL_RNA_646 TS047_RoL_RNA_651 TS047_RoL_RNA_652 TS047_RoL_RNA_655 TS047_RoL_RNA_676 TS047_RoL_RNA_677 TS047_RoL_RNA_678 TS047_RoL_RNA_679 TS047_RoL_RNA_683 TS047_RoL_RNA_707 TS047_RoL_RNA_708 TS047_RoL_RNA_710 TS047_RoL_RNA_79 TS047_RoL_RNA_85

cat("Available assays:", assayNames(se), "\n")

## Available assays: counts abundance

cat("Column data (metadata):", names(colData(se)), "\n")

## Column data (metadata): sample fastq_1 fastq_2 strandedness

Data Preprocessing

# Clean sample names by removing "TS047_RoL_RNA_" prefix
clean_names <- gsub("TS047_RoL_RNA_", "", colnames(se))

# Update column names in the SummarizedExperiment object
colnames(se) <- clean_names

# Create a data frame for the metadata columns we want to add
metadata_cols <- c("Time", "Antibiotics", "Temperature", 
                  "Pathogen", "History", "Length.mm", "Weight.g", "Total.Worm.Count")

# Initialize a data frame to store the metadata
sample_metadata <- data.frame(
  Sample_Name = clean_names,
  stringsAsFactors = FALSE
)
rownames(sample_metadata) <- clean_names

# Add each metadata column
for (col in metadata_cols) {
  values <- rep(NA, length(clean_names))
  
  for (i in seq_along(clean_names)) {
    matches <- which(Metadata$Sample_Name == clean_names[i])
    
    if (length(matches) == 1) {
      values[i] <- Metadata[matches, col, drop = TRUE]
    }
  }
  
  sample_metadata[[col]] <- values
}

# Add the metadata to the SummarizedExperiment object
colData(se) <- DataFrame(sample_metadata)

cat("Metadata added to SummarizedExperiment object\n")

## Metadata added to SummarizedExperiment object

Create DESeq2 Dataset

# Remove samples with NA values in treatment columns
se_filtered <- se[, !is.na(colData(se)$Antibiotics) & !is.na(colData(se)$Temperature) & !is.na(colData(se)$Pathogen)]
cat("Number of samples after removing NAs:", ncol(se_filtered), "\n")

## Number of samples after removing NAs: 72

# Create treatment combinations using DESeq2-compatible format
colData(se_filtered)$Treatment <- paste0(
  ifelse(colData(se_filtered)$Antibiotics == 1, "A_plus", "A_minus"), "_",
  ifelse(colData(se_filtered)$Temperature == 1, "T_plus", "T_minus"), "_",
  ifelse(colData(se_filtered)$Pathogen == 1, "P_plus", "P_minus")
)

# Create plotting-friendly treatment names
colData(se_filtered)$Treatment_Plot <- paste0(
  ifelse(colData(se_filtered)$Antibiotics == 1, "A+", "A-"), " ",
  ifelse(colData(se_filtered)$Temperature == 1, "T+", "T-"), " ",
  ifelse(colData(se_filtered)$Pathogen == 1, "P+", "P-")
)

# Convert treatment to factor with specified order (DESeq2 names)
colData(se_filtered)$Treatment <- factor(colData(se_filtered)$Treatment, levels = treatment_order_deseq)

# Convert plotting treatment to factor with specified order
colData(se_filtered)$Treatment_Plot <- factor(colData(se_filtered)$Treatment_Plot, levels = treatment_order_plot)

# Print the distribution of treatment combinations
cat("\nTreatment combination distribution:\n")

## 
## Treatment combination distribution:

print(table(colData(se_filtered)$Treatment))

## 
## A_minus_T_minus_P_minus  A_minus_T_minus_P_plus  A_plus_T_minus_P_minus 
##                       6                      12                       6 
##   A_plus_T_minus_P_plus  A_minus_T_plus_P_minus   A_minus_T_plus_P_plus 
##                      12                       6                      12 
##   A_plus_T_plus_P_minus    A_plus_T_plus_P_plus 
##                       6                      12

print(table(colData(se_filtered)$Treatment_Plot))

## 
## A- T- P- A- T- P+ A+ T- P- A+ T- P+ A- T+ P- A- T+ P+ A+ T+ P- A+ T+ P+ 
##        6       12        6       12        6       12        6       12

# Clean other factor levels
colData(se_filtered)$Time <- factor(colData(se_filtered)$Time)
colData(se_filtered)$Antibiotics <- factor(colData(se_filtered)$Antibiotics)
colData(se_filtered)$Temperature <- factor(colData(se_filtered)$Temperature)
colData(se_filtered)$Pathogen <- factor(colData(se_filtered)$Pathogen)

# Create initial design formula
design_formula <- ~ Treatment

# We need to use raw counts for DESeq2, not length-scaled counts
if (!"counts" %in% assayNames(se_filtered)) {
  # Read in the raw counts file
  raw_counts <- read.table(here::here("Data", "Transcriptomics", "Results", "rnaseq", "051525", "salmon", "salmon.merged.gene_counts_length_scaled.tsv"), 
                          header = TRUE, row.names = 1)
  
  # Convert to matrix and ensure it's numeric
  counts_matrix <- as.matrix(raw_counts)
  mode(counts_matrix) <- "numeric"
  
  # Round to nearest integer (DESeq2 requires integer counts)
  counts_matrix <- round(counts_matrix)
  
  # Ensure the order of samples matches the SummarizedExperiment
  counts_matrix <- counts_matrix[, colnames(se_filtered)]
} else {
  # If counts are available, use them
  counts_matrix <- as.matrix(assays(se_filtered)$counts)
  mode(counts_matrix) <- "numeric"
  counts_matrix <- round(counts_matrix)
}

# Create initial DESeq2 dataset
dds_initial <- DESeqDataSetFromMatrix(
  countData = counts_matrix,
  colData = colData(se_filtered),
  design = design_formula
)

## converting counts to integer mode

# Pre-filtering: remove low abundance transcripts
cpm <- counts(dds_initial) / colSums(counts(dds_initial)) * 1e6

# Keep genes that have at least 2 CPM in at least 3 samples
keep <- rowSums(cpm >= 2) >= 3

dds_filtered <- dds_initial[keep,]

# Print filtering statistics
cat("Number of genes before filtering:", nrow(counts_matrix), "\n")

## Number of genes before filtering: 39447

cat("Number of genes after filtering:", nrow(dds_filtered), "\n")

## Number of genes after filtering: 20672

cat("Number of genes removed:", nrow(counts_matrix) - nrow(dds_filtered), "\n")

## Number of genes removed: 18775

cat("Percentage of genes kept:", round(nrow(dds_filtered)/nrow(counts_matrix)*100, 2), "%\n")

## Percentage of genes kept: 52.4 %

# Transform counts for visualization
rld <- vst(dds_filtered, blind=FALSE)

# Estimate size factors
dds_filtered <- estimateSizeFactors(dds_filtered)

# Get normalized counts
normalized_counts <- counts(dds_filtered, normalized=TRUE)

# After dds_filtered is created and before any downstream analyses:

# List of genes to exclude
genes_to_exclude <- c(
  "zgc:112970",
  "lectin",
  "zp2.6",
  "qdprb.4",
  "bncr",
  "avd",
  "tdgf1",
  "adad2",
  "LOC137488852",
  "ppt2a.1",
  "snapc1a",
  "exd1",
  "LOC141381090",
  "LOC141381093",
  "LOC141378343",
  "LOC101884348"
)

# Remove these genes from dds_filtered, rld, and normalized_counts
keep_genes <- !(rownames(dds_filtered) %in% genes_to_exclude)
dds_filtered <- dds_filtered[keep_genes, ]

# Re-transform counts for visualization after filtering
rld <- vst(dds_filtered, blind=FALSE)

# Estimate size factors again after filtering
dds_filtered <- estimateSizeFactors(dds_filtered)

# Get normalized counts after filtering
normalized_counts <- counts(dds_filtered, normalized=TRUE)

# Print exclusion summary
cat("Number of genes excluded:", length(genes_to_exclude), "\n")

## Number of genes excluded: 16

cat("Number of genes remaining after exclusion:", nrow(dds_filtered), "\n")

## Number of genes remaining after exclusion: 20657

Set up parallel processing

# Enhanced parallel processing setup
# Detect number of available cores
available_cores <- parallel::detectCores()
cat("Available cores:", available_cores, "\n")

## Available cores: 10

# Use 75% of available cores to avoid overwhelming the system
# Leave some cores free for other processes
num_cores <- max(1, floor(available_cores * 0.75))
cat("Using cores for parallel processing:", num_cores, "\n")

## Using cores for parallel processing: 7

# Create optimized BiocParallel parameters
bp_param <- MulticoreParam(
  workers = num_cores,
  tasks = num_cores * 2,  # Increase task granularity
  progressbar = TRUE      # Enable progress reporting
)

# Register the optimized parameters
register(bp_param)

cat("Parallel processing configured with", num_cores, "cores\n")

## Parallel processing configured with 7 cores

cat("Task granularity:", bp_param$tasks, "\n")

## Task granularity: 14

cat("Progress bar enabled:", bp_param$progressbar, "\n")

## Progress bar enabled: TRUE

Define Analysis Comparisons

# Define the specific comparisons for our analysis
comparisons <- list(
  
  # 1. All Treatment Groups (comprehensive analysis)
  "All_Treatments" = list(
    description = "Comprehensive analysis across all treatment combinations",
    design = ~ Treatment,
    contrasts = list(
      "A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_minus_T_minus_P_plus", "A_minus_T_minus_P_minus"),
      "A_plus_T_minus_P_minus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_plus_T_minus_P_minus", "A_minus_T_minus_P_minus"),
      "A_plus_T_minus_P_plus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_plus_T_minus_P_plus", "A_minus_T_minus_P_minus"),
      "A_minus_T_plus_P_minus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_minus_T_plus_P_minus", "A_minus_T_minus_P_minus"),
      "A_minus_T_plus_P_plus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_minus_T_plus_P_plus", "A_minus_T_minus_P_minus"),
      "A_plus_T_plus_P_minus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_plus_T_plus_P_minus", "A_minus_T_minus_P_minus"),
      "A_plus_T_plus_P_plus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_plus_T_plus_P_plus", "A_minus_T_minus_P_minus")
    )
  ),
  
  # 2. Parasite Effect (A_minus_T_minus_P_minus vs A_minus_T_minus_P_plus)
  "Parasite_Effect" = list(
    description = "Baseline parasite effect: A_minus_T_minus_P_minus vs A_minus_T_minus_P_plus",
    design = ~ Treatment,
    contrasts = list(
      "A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_minus_T_minus_P_plus", "A_minus_T_minus_P_minus")
    )
  ),
  
  # 3. Historical Contingency (A_minus_T_minus_P_plus vs A_plus_T_minus_P_plus, A_minus_T_plus_P_plus, A_plus_T_plus_P_plus)
  "Historical_Contingency" = list(
    description = "How prior stressors affect parasite response: A_minus_T_minus_P_plus vs (A_plus_T_minus_P_plus, A_minus_T_plus_P_plus, A_plus_T_plus_P_plus)",
    design = ~ Treatment,
    contrasts = list(
      "A_plus_T_minus_P_plus_vs_A_minus_T_minus_P_plus" = c("Treatment", "A_plus_T_minus_P_plus", "A_minus_T_minus_P_plus"),
      "A_minus_T_plus_P_plus_vs_A_minus_T_minus_P_plus" = c("Treatment", "A_minus_T_plus_P_plus", "A_minus_T_minus_P_plus"),
      "A_plus_T_plus_P_plus_vs_A_minus_T_minus_P_plus" = c("Treatment", "A_plus_T_plus_P_plus", "A_minus_T_minus_P_plus")
    )
  ),
  
  # 4. Recovery Analysis (A_minus_T_minus_P_minus vs A_plus_T_minus_P_minus, A_minus_T_plus_P_minus, A_plus_T_plus_P_minus)
  "Recovery_Analysis" = list(
    description = "How prior stressors affect recovery: A_minus_T_minus_P_minus vs (A_plus_T_minus_P_minus, A_minus_T_plus_P_minus, A_plus_T_plus_P_minus)",
    design = ~ Treatment,
    contrasts = list(
      "A_plus_T_minus_P_minus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_plus_T_minus_P_minus", "A_minus_T_minus_P_minus"),
      "A_minus_T_plus_P_minus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_minus_T_plus_P_minus", "A_minus_T_minus_P_minus"),
      "A_plus_T_plus_P_minus_vs_A_minus_T_minus_P_minus" = c("Treatment", "A_plus_T_plus_P_minus", "A_minus_T_minus_P_minus")
    )
  )
)

Code

(hide)

Quality Control

Sample Distance Heatmap

Principal Component Analysis

## using ntop=500 top features by variance

Manage Saved Results

# Check which analyses have saved results
cat("=== Checking for saved results ===\n")

## === Checking for saved results ===

for (analysis_name in names(comparisons)) {
  if (results_exist(analysis_name)) {
    cat("✓", analysis_name, "- results found\n")
  } else {
    cat("✗", analysis_name, "- no saved results\n")
  }
}

## ✓ All_Treatments - results found
## ✓ Parasite_Effect - results found
## ✓ Historical_Contingency - results found
## ✓ Recovery_Analysis - results found

# Set this to TRUE if you want to recompute all analyses
# Set to FALSE to use existing results when available
force_recompute_all <- FALSE

# You can also set individual analyses to recompute
force_recompute_individual <- list(
  "All_Treatments" = FALSE,
  "Parasite_Effect" = FALSE,
  "Historical_Contingency" = FALSE,
  "Recovery_Analysis" = FALSE
)

cat("\nForce recompute all:", force_recompute_all, "\n")

## 
## Force recompute all: FALSE

cat("Individual force recompute settings:\n")

## Individual force recompute settings:

for (analysis_name in names(force_recompute_individual)) {
  cat("  ", analysis_name, ":", force_recompute_individual[[analysis_name]], "\n")
}

##    All_Treatments : FALSE 
##    Parasite_Effect : FALSE 
##    Historical_Contingency : FALSE 
##    Recovery_Analysis : FALSE

# Uncomment the lines below if you want to clean up saved results
# cleanup_saved_results()  # Remove all saved results
# cleanup_saved_results(c("All_Treatments", "Parasite_Effect"))  # Remove specific analyses

cat("\n=== Usage Instructions ===\n")

## 
## === Usage Instructions ===

cat("1. Set force_recompute_all = TRUE to recompute all analyses\n")

## 1. Set force_recompute_all = TRUE to recompute all analyses

cat("2. Set individual force_recompute_individual values to TRUE to recompute specific analyses\n")

## 2. Set individual force_recompute_individual values to TRUE to recompute specific analyses

cat("3. Use cleanup_saved_results() to remove all saved results\n")

## 3. Use cleanup_saved_results() to remove all saved results

cat("4. Use cleanup_saved_results(c('analysis1', 'analysis2')) to remove specific analyses\n")

## 4. Use cleanup_saved_results(c('analysis1', 'analysis2')) to remove specific analyses

Run DESeq2 Analyses

# Monitor initial performance
cat("=== Starting DESeq2 Analysis ===\n")

## === Starting DESeq2 Analysis ===

monitor_performance()

## === Performance Monitor ===
## Memory usage:
##            used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells  7248456 387.2   11112538 593.5         NA 11112538 593.5
## Vcells 29991555 228.9   54056739 412.5      32768 54056739 412.5
## CPU info:
## CPU usage: 37.15% user, 15.20% sys, 47.63% idle  
## Active parallel workers: 7 
## ========================

# Optimize memory before starting
optimize_memory()

## Optimizing memory usage...
## Memory optimization complete

# Run all analyses
all_results <- list()

for (analysis_name in names(comparisons)) {
  cat("\n--- Starting analysis:", analysis_name, "---\n")
  
  # Monitor performance before each analysis
  monitor_performance()
  
  # Run the analysis
  all_results[[analysis_name]] <- run_analysis_set(analysis_name, comparisons[[analysis_name]], force_recompute_all || force_recompute_individual[[analysis_name]])
  
  # Optimize memory after each analysis
  optimize_memory()
  
  cat("--- Completed analysis:", analysis_name, "---\n")
}

## 
## --- Starting analysis: All_Treatments ---
## === Performance Monitor ===
## Memory usage:
##            used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells  7248618 387.2   11112538 593.5         NA 11112538 593.5
## Vcells 29992406 228.9   54056739 412.5      32768 54056739 412.5
## CPU info:
## CPU usage: 35.8% user, 14.66% sys, 50.24% idle  
## Active parallel workers: 7 
## ========================
## 
## === Running All_Treatments ===
## Description: Comprehensive analysis across all treatment combinations 
## Loading existing results for All_Treatments 
## Loading existing results for: All_Treatments 
## Optimizing memory usage...
## Memory optimization complete
## --- Completed analysis: All_Treatments ---
## 
## --- Starting analysis: Parasite_Effect ---
## === Performance Monitor ===
## Memory usage:
##            used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells  7249477 387.2   11112538 593.5         NA 11112538 593.5
## Vcells 31010273 236.6   54056739 412.5      32768 54056739 412.5
## CPU info:
## CPU usage: 35.67% user, 12.16% sys, 52.16% idle  
## Active parallel workers: 7 
## ========================
## 
## === Running Parasite_Effect ===
## Description: Baseline parasite effect: A_minus_T_minus_P_minus vs A_minus_T_minus_P_plus 
## Loading existing results for Parasite_Effect 
## Loading existing results for: Parasite_Effect 
## Optimizing memory usage...
## Memory optimization complete
## --- Completed analysis: Parasite_Effect ---
## 
## --- Starting analysis: Historical_Contingency ---
## === Performance Monitor ===
## Memory usage:
##            used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells  7249616 387.2   11112538 593.5         NA 11112538 593.5
## Vcells 31155745 237.7   54056739 412.5      32768 54056739 412.5
## CPU info:
## CPU usage: 40.50% user, 12.3% sys, 47.45% idle  
## Active parallel workers: 7 
## ========================
## 
## === Running Historical_Contingency ===
## Description: How prior stressors affect parasite response: A_minus_T_minus_P_plus vs (A_plus_T_minus_P_plus, A_minus_T_plus_P_plus, A_plus_T_plus_P_plus) 
## Loading existing results for Historical_Contingency 
## Loading existing results for: Historical_Contingency 
## Optimizing memory usage...
## Memory optimization complete
## --- Completed analysis: Historical_Contingency ---
## 
## --- Starting analysis: Recovery_Analysis ---
## === Performance Monitor ===
## Memory usage:
##            used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells  7249841 387.2   11112538 593.5         NA 11112538 593.5
## Vcells 31591152 241.1   54056739 412.5      32768 54056739 412.5
## CPU info:
## CPU usage: 37.11% user, 14.51% sys, 48.36% idle  
## Active parallel workers: 7 
## ========================
## 
## === Running Recovery_Analysis ===
## Description: How prior stressors affect recovery: A_minus_T_minus_P_minus vs (A_plus_T_minus_P_minus, A_minus_T_plus_P_minus, A_plus_T_plus_P_minus) 
## Loading existing results for Recovery_Analysis 
## Loading existing results for: Recovery_Analysis 
## Optimizing memory usage...
## Memory optimization complete
## --- Completed analysis: Recovery_Analysis ---

# Final performance check
cat("\n=== Final Performance Check ===\n")

## 
## === Final Performance Check ===

monitor_performance()

## === Performance Monitor ===
## Memory usage:
##            used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells  7250003 387.2   11112538 593.5         NA 11112538 593.5
## Vcells 32026015 244.4   54056739 412.5      32768 54056739 412.5
## CPU info:
## CPU usage: 35.85% user, 13.92% sys, 50.21% idle  
## Active parallel workers: 7 
## ========================

# Print summary of results
cat("\n=== Analysis Summary ===\n")

## 
## === Analysis Summary ===

for (analysis_name in names(all_results)) {
  cat("\n", analysis_name, ":\n")
  for (contrast_name in names(all_results[[analysis_name]])) {
    res <- all_results[[analysis_name]][[contrast_name]]
    cat("  ", contrast_name, ":\n")
    cat("    Significant genes (padj < 0.05):", sum(res$padj < 0.05, na.rm=TRUE), "\n")
    cat("    Upregulated genes:", sum(res$padj < 0.05 & res$log2FoldChange > 0, na.rm=TRUE), "\n")
    cat("    Downregulated genes:", sum(res$padj < 0.05 & res$log2FoldChange < 0, na.rm=TRUE), "\n")
  }
}

## 
##  All_Treatments :
##    A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 3822 
##     Upregulated genes: 2012 
##     Downregulated genes: 1810 
##    A_plus_T_minus_P_minus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 42 
##     Upregulated genes: 12 
##     Downregulated genes: 30 
##    A_plus_T_minus_P_plus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 6411 
##     Upregulated genes: 3426 
##     Downregulated genes: 2985 
##    A_minus_T_plus_P_minus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 66 
##     Upregulated genes: 13 
##     Downregulated genes: 53 
##    A_minus_T_plus_P_plus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 6144 
##     Upregulated genes: 3213 
##     Downregulated genes: 2931 
##    A_plus_T_plus_P_minus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 123 
##     Upregulated genes: 22 
##     Downregulated genes: 101 
##    A_plus_T_plus_P_plus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 4930 
##     Upregulated genes: 2574 
##     Downregulated genes: 2356 
## 
##  Parasite_Effect :
##    A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 3822 
##     Upregulated genes: 2012 
##     Downregulated genes: 1810 
## 
##  Historical_Contingency :
##    A_plus_T_minus_P_plus_vs_A_minus_T_minus_P_plus :
##     Significant genes (padj < 0.05): 55 
##     Upregulated genes: 33 
##     Downregulated genes: 22 
##    A_minus_T_plus_P_plus_vs_A_minus_T_minus_P_plus :
##     Significant genes (padj < 0.05): 132 
##     Upregulated genes: 95 
##     Downregulated genes: 37 
##    A_plus_T_plus_P_plus_vs_A_minus_T_minus_P_plus :
##     Significant genes (padj < 0.05): 22 
##     Upregulated genes: 6 
##     Downregulated genes: 16 
## 
##  Recovery_Analysis :
##    A_plus_T_minus_P_minus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 42 
##     Upregulated genes: 12 
##     Downregulated genes: 30 
##    A_minus_T_plus_P_minus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 66 
##     Upregulated genes: 13 
##     Downregulated genes: 53 
##    A_plus_T_plus_P_minus_vs_A_minus_T_minus_P_minus :
##     Significant genes (padj < 0.05): 123 
##     Upregulated genes: 22 
##     Downregulated genes: 101

Results Display

Heatmap

## Number of genes excluded: 16 
## Number of genes remaining in top_genes_list: 97

00) All Treatment Groups Analysis

Volcano Plots

Summary Tables

summary_table

Treatment Comparison	Total Genes	Significant Genes	Upregulated	Downregulated
Summary of Differential Gene Expression
All Treatment Groups Analysis
A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus	20,672	3,822	2,012	1,810
A_plus_T_minus_P_minus_vs_A_minus_T_minus_P_minus	20,672	42	12	30
A_plus_T_minus_P_plus_vs_A_minus_T_minus_P_minus	20,672	6,411	3,426	2,985
A_minus_T_plus_P_minus_vs_A_minus_T_minus_P_minus	20,672	66	13	53
A_minus_T_plus_P_plus_vs_A_minus_T_minus_P_minus	20,672	6,144	3,213	2,931
A_plus_T_plus_P_minus_vs_A_minus_T_minus_P_minus	20,672	123	22	101
A_plus_T_plus_P_plus_vs_A_minus_T_minus_P_minus	20,672	4,930	2,574	2,356

Top Genes Tables

## ### Top 10 Upregulated Genes: A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Downregulated Genes: A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Upregulated Genes: A_plus_T_minus_P_minus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Downregulated Genes: A_plus_T_minus_P_minus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Upregulated Genes: A_plus_T_minus_P_plus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Downregulated Genes: A_plus_T_minus_P_plus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Upregulated Genes: A_minus_T_plus_P_minus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Downregulated Genes: A_minus_T_plus_P_minus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Upregulated Genes: A_minus_T_plus_P_plus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Downregulated Genes: A_minus_T_plus_P_plus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Upregulated Genes: A_plus_T_plus_P_minus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Downregulated Genes: A_plus_T_plus_P_minus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Upregulated Genes: A_plus_T_plus_P_plus_vs_A_minus_T_minus_P_minus 
## ### Top 10 Downregulated Genes: A_plus_T_plus_P_plus_vs_A_minus_T_minus_P_minus

01) Parasite Effect Analysis

Volcano Plot

Summary Table

summary_table

Treatment Comparison	Total Genes	Significant Genes	Upregulated	Downregulated
Summary of Differential Gene Expression
Parasite Effect Analysis (A- T- P- vs A- T- P+)
A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus	20,672	3,822	2,012	1,810

Top Genes Table

## ### Top 20 Significant Genes: A_minus_T_minus_P_plus_vs_A_minus_T_minus_P_minus

02) Historical Contingency Analysis

Volcano Plots

Summary Table

summary_table

Treatment Comparison	Upregulated Genes	Downregulated Genes	Total Significant
Historical Contingency Summary
Effect of prior stressors on parasite response
ABX+ vs Control	33	22	55
TEMP+ vs Control	95	37	132
ABX+TEMP+ vs Control	6	16	22

Top Genes Tables

## ### Top 15 Significant Genes: A_plus_T_minus_P_plus_vs_A_minus_T_minus_P_plus 
## ### Top 15 Significant Genes: A_minus_T_plus_P_plus_vs_A_minus_T_minus_P_plus 
## ### Top 15 Significant Genes: A_plus_T_plus_P_plus_vs_A_minus_T_minus_P_plus

03) Recovery Analysis

Volcano Plots

Summary Table

summary_table

Treatment Comparison	Upregulated Genes	Downregulated Genes	Total Significant
Recovery Analysis Summary
Effect of prior stressors on recovery
ABX+ vs Control	12	30	42
TEMP+ vs Control	13	53	66
ABX+TEMP+ vs Control	22	101	123

Top Genes Tables

## ### Top 15 Significant Genes: A_plus_T_minus_P_minus_vs_A_minus_T_minus_P_minus 
## ### Top 15 Significant Genes: A_minus_T_plus_P_minus_vs_A_minus_T_minus_P_minus 
## ### Top 15 Significant Genes: A_plus_T_plus_P_minus_vs_A_minus_T_minus_P_minus

Summary Plots

All Treatment Groups Analysis

summary_table

Treatment Comparison	Upregulated Genes	Downregulated Genes	Total Significant
Summary Statistics - All Treatment Groups
Number of significantly differentially expressed genes
A- T- P+ vs. A- T- P-	2,012	1,810	3,822
A+ T- P- vs. A- T- P-	12	30	42
A+ T- P+ vs. A- T- P-	3,426	2,985	6,411
A- T+ P- vs. A- T- P-	13	53	66
A- T+ P+ vs. A- T- P-	3,213	2,931	6,144
A+ T+ P- vs. A- T- P-	22	101	123
A+ T+ P+ vs. A- T- P-	2,574	2,356	4,930

Parasite Effect Analysis

summary_table

Treatment Comparison	Upregulated Genes	Downregulated Genes	Total Significant	PlotName
Summary Statistics - Parasite Effect
Number of significantly differentially expressed genes
A- T- P+ vs. A- T- P-	2,012	1,810	3,822	A- T- P+ vs. A- T- P-

Historical Contingency Analysis

summary_table

Treatment Comparison	Upregulated Genes	Downregulated Genes	Total Significant
Historical Contingency Summary
Effect of prior stressors on parasite response
A+ T- P+ vs. A- T- P-	33	22	55
A- T+ P+ vs. A- T- P-	95	37	132
A+ T+ P+ vs. A- T- P-	6	16	22

Recovery Analysis

summary_table

Treatment Comparison	Upregulated Genes	Downregulated Genes	Total Significant
Recovery Analysis Summary
Effect of prior stressors on recovery
A+ T- P- vs. A- T- P-	12	30	42
A- T+ P- vs. A- T- P-	13	53	66
A+ T+ P- vs. A- T- P-	22	101	123

Combined Summary Across All Analyses

comprehensive_table

Analysis Type	Treatment Combination	Upregulated Genes	Downregulated Genes	Total Significant
Comprehensive Differential Gene Expression Summary
All analyses and comparisons
All_Treatments	A- T- P+ vs A- T- P-	2,012	1,810	3,822
All_Treatments	ABX vs A- T- P-	12	30	42
All_Treatments	A+ T- P+ vs A- T- P-	3,426	2,985	6,411
All_Treatments	A- T+ P- vs A- T- P-	13	53	66
All_Treatments	A- T+ P+ vs A- T- P-	3,213	2,931	6,144
All_Treatments	A+ T+ P- vs A- T- P-	22	101	123
All_Treatments	A+ T+ P+ vs A- T- P-	2,574	2,356	4,930
Parasite_Effect	A- T- P+ vs A- T- P-	2,012	1,810	3,822
Historical_Contingency	A+ T- P+ vs A- T- P+	33	22	55
Historical_Contingency	A- T+ P+ vs A- T- P+	95	37	132
Historical_Contingency	A+ T+ P+ vs A- T- P+	6	16	22
Recovery_Analysis	ABX vs A- T- P-	12	30	42
Recovery_Analysis	A- T+ P- vs A- T- P-	13	53	66
Recovery_Analysis	A+ T+ P- vs A- T- P-	22	101	123

Export Results

Session Info

sessionInfo()

## R version 4.3.3 (2024-02-29)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS 15.3
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] coin_1.4-3                  rstatix_0.7.2              
##  [3] multcomp_1.4-25             TH.data_1.1-2              
##  [5] MASS_7.3-60.0.1             survival_3.6-4             
##  [7] mvtnorm_1.2-5               emmeans_1.10.1             
##  [9] car_3.1-2                   carData_3.0-5              
## [11] broom_1.0.6                 gt_0.10.1                  
## [13] BiocParallel_1.36.0         gridExtra_2.3              
## [15] VennDiagram_1.7.3           futile.logger_1.4.3        
## [17] lubridate_1.9.3             forcats_1.0.0              
## [19] stringr_1.5.1               dplyr_1.1.4                
## [21] purrr_1.0.2                 readr_2.1.5                
## [23] tidyr_1.3.1                 tibble_3.2.1               
## [25] tidyverse_2.0.0             ggrepel_0.9.5              
## [27] RColorBrewer_1.1-3          pheatmap_1.0.12            
## [29] ggplot2_3.5.1               DESeq2_1.42.1              
## [31] SummarizedExperiment_1.32.0 Biobase_2.62.0             
## [33] MatrixGenerics_1.14.0       matrixStats_1.3.0          
## [35] GenomicRanges_1.54.1        GenomeInfoDb_1.38.8        
## [37] IRanges_2.36.0              S4Vectors_0.40.2           
## [39] BiocGenerics_0.48.1        
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-7            formatR_1.14            readxl_1.4.3           
##  [4] sandwich_3.1-0          rlang_1.1.4             magrittr_2.0.3         
##  [7] compiler_4.3.3          png_0.1-8               vctrs_0.6.5            
## [10] pkgconfig_2.0.3         crayon_1.5.2            fastmap_1.2.0          
## [13] backports_1.5.0         XVector_0.42.0          labeling_0.4.3         
## [16] utf8_1.2.4              rmarkdown_2.26          tzdb_0.4.0             
## [19] modeltools_0.2-23       xfun_0.43               zlibbioc_1.48.2        
## [22] cachem_1.1.0            jsonlite_1.8.8          highr_0.10             
## [25] DelayedArray_0.28.0     parallel_4.3.3          R6_2.5.1               
## [28] bslib_0.7.0             stringi_1.8.4           cellranger_1.1.0       
## [31] jquerylib_0.1.4         estimability_1.5        Rcpp_1.0.12            
## [34] knitr_1.45              zoo_1.8-12              Matrix_1.6-5           
## [37] splines_4.3.3           timechange_0.3.0        tidyselect_1.2.1       
## [40] rstudioapi_0.16.0       abind_1.4-5             yaml_2.3.8             
## [43] codetools_0.2-20        lattice_0.22-6          withr_3.0.0            
## [46] coda_0.19-4.1           evaluate_0.23           lambda.r_1.2.4         
## [49] xml2_1.3.6              pillar_1.9.0            generics_0.1.3         
## [52] rprojroot_2.0.4         RCurl_1.98-1.14         hms_1.1.3              
## [55] munsell_0.5.1           scales_1.3.0            xtable_1.8-4           
## [58] glue_1.7.0              tools_4.3.3             locfit_1.5-9.9         
## [61] libcoin_1.0-10          colorspace_2.1-0        GenomeInfoDbData_1.2.11
## [64] cli_3.6.2               futile.options_1.0.1    fansi_1.0.6            
## [67] S4Arrays_1.2.1          gtable_0.3.5            sass_0.4.9             
## [70] digest_0.6.35           SparseArray_1.2.4       farver_2.1.2           
## [73] htmltools_0.5.8.1       lifecycle_1.0.4         here_1.0.1

Historical Contingency of Host-Microbiome Response to Perturbation

Differential Gene Expression Results

Michael Sieler

28 June, 2025

Overview

Setup

(hide)

Libraries

Plotting

Functions

Import Data

Data Preprocessing

Create DESeq2 Dataset

Set up parallel processing

Define Analysis Comparisons

Code

(hide)

Quality Control

Sample Distance Heatmap

Principal Component Analysis

Manage Saved Results

Run DESeq2 Analyses

Results Display

Heatmap

00) All Treatment Groups Analysis

Volcano Plots

Summary Tables

Top Genes Tables

01) Parasite Effect Analysis

Volcano Plot

Summary Table

Top Genes Table

02) Historical Contingency Analysis

Volcano Plots

Summary Table

Top Genes Tables

03) Recovery Analysis

Volcano Plots

Summary Table

Top Genes Tables

Summary Plots

All Treatment Groups Analysis

Parasite Effect Analysis

Historical Contingency Analysis

Recovery Analysis

Combined Summary Across All Analyses

Export Results

Session Info